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Abstract 

The implementation of the kinetic diffuse boundary condition with the characteristic streaming- 
collision mechanism is studied for the high-order lattice Boltzmann (LB) models. The obtained 
formulation is also tested and validated numerically for three high-order LB models for both isother- 
mal and thermal Couette flows. The streaming-collision mechanism ensures that high-order LB 
models can retain particle feature while go beyond the Navier-Stokes hydrodynamics. 
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I. INTRODUCTION 



High-order models have recently attracted considerable interests in the lattice Boltzmann 
(LB) community. For these models, high-order terms in the expanded distribution function 
are retained, thus multi-speed lattices have to be used. By doing so, we have a few ben- 
efits, e.g. consistent description of thermal flows, the Galilean invariance of the transport 
coefficients, improved model capability for compressible and rarefied flowsp~H3]. Meanwhile, 
high-order models can still preserve the simplicity of the standard LB model. High-order LB 
models are tererfore often applied to thermal flows, compressible flows, and rarefied flows 

ESI El- 

Due to their kinetic origin, high-order LB models have shown to be able to approach 
the Boltzmann model equation such as the Bhatnagar-Gross-Krook (BGK) equation by 
increasing the expansion and quadrature order pQ. Numerically, it is found that high-order 
models are capable of describing rarefaction effects up to the early transitional regime with 
slightly increased discrete velocities. For example, the velocity-slip can be captured by a 
D2Q16 model for the Knudsen number up to 0.5 where DmQn refers to the standard m 
dimenional and n discrete velocities. A recent study on thermal flows shows that thermal 
rarefaction effects can also be captured by using a moderate number of discrete velocities 
with an affordable computational cost (5]. In addition to capturing the non-equilibrium 
effects, high-order models are also useful for other applications. For instance, by using a 
D2Q25 model, Daniel et at [2] have shown that high Weber number can be achievable for 
droplet collision simulations. 

Implementation of boundary conditions is crucial to application of high-order LB models. 
The kinetic diffuse-reflection boundary condition has shown to be able to predict velocity- 
slip and temperature-jump at the solid boundary. Moreover, the positivity of distribution 
function can always be maintained (provided that the distribution function from the bulk 
is positive), which is key to numerical stability. However, due to the multi-speed lattice 
of high-order models, the formulation of the kinetic diffuse-reflection boundary condition 
with the characteristic "streaming and collision" mechanism is yet to be developed. The 
successful implementations so far for high-order models are based on various finite difference 
scheme [8], where the highly desirable "streaming and collision" mechanism disappears. 
Because of this "streaming and collision" mechanism, LB method is often regarded a particle 
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method. The purpose of this work is to formulate the kinetic diffuse reflection boundary 
condition to retain the "streaming and collision" feature. 



II. BRIEF DESCRIPTION OF HIGH ORDER LATTICE BOLTZMANN MODEL 

The Boltzmann-BGK equation can be dicretized using a systematic procedure (see [TJ El 
[TO] for the details) to derive the LB governing equation, which can be written as 
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where / denotes the single-particle distribution function evaluated at a discrete velocity 
c a , f eq is the truncated Maxwellian distribution, while r is the mean relaxation time. For 
convenience, a non-dimensional system 
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can be introduced in Eq. ([!]), where the symbols with hat represent the dimensional quanti- 
ties. The common notations are used to represent physical quantities, i.e. p denotes density; 
/i, dynamic viscosity; u, velocity; p, pressure; T, temperature; a, stress and q, heat flux. L 
is the characteristic length of the system. The symbols with subscript are the correspond- 
ing reference quantities. With this non-dimensional system, the equation of state becomes 
p = pT. The mean relaxation time can be written explicitly as r = po^RTop/ (p Lp), 
which is related to the viscosity and pressure. Meanwhile, the macroscopic quantities can 
be obtained as 
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where D is the space dimension number and d is the total discrete velocity number and 
the angle brackets < ■ ■ ■ > indicates the trace-free part of the tensor. The discrete veloc- 
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ity c a and its weights w a may be determined through several ways, e.g., [H [TTJ [12] list 
explicitly various orders of discrete velocity sets. For convenience, we use the notation 
£ = {c a ,w a },a = l..d to represent the discrete velocity set. With an appropriate £, the 
explicit form of the truncated Maxwellian distribution can be given as 
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in which the forth order expansion is used. If an isothermal flow is concerned, the tempera- 
ture T should be set to 1. Moreover, to simulate incompressible flows, it is common to use 
only the second order terms of g^. 

With Eq.(|T| and Q, the final issue is to choose an appropriate numerical scheme. Fol- 
lowing the spirit of "streaming and collision" mechanism [13J, an implicit scheme 
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dt 



2t(x + c a dt, t + dt) 
can be constructed. If we introduce a new variable 

fa = fa~\~ ^2/j-^ a ~ /a 9 ) 

to eliminate the implicitness, we will get the evolution equation for/ as 



[fa q (x + c a dt, t + dt) - f a (x + Cadt, t + dt)} 



f a (x + C a dt,t + dt) - f a (x,t) 
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t(x, t) + 0.5dt 

The advantage of Eq.(|5]) is that, if the discrete velocities are tied to discretization of the 
space and time by choosing £ with integer value, the evolution of f a can be accomplished 
in a way similar to a "particle", which makes the LB method simple but still flexible. 
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With the variable f a , conservative quantities like density can still be obtained by using 
Eq.(j3]) without changing form but some conversions are needed for shear stress and heat 
flux (cf. [13]). In addition, the mean relaxation time r may be related to the local gas 
temperature for thermal problems. 

III. KINETIC DIFFUSE-REFLECTION-TYPE BOUNDARY CONDITION 

The essential idea of the kinetic diffuse-reflection boundary condition is that an outgoing 
particle completely forgets its history and its velocity is re-normalized by the Maxwellian 
distribution. To implement this diffuse-reflection principle for the high-order LB model, 
we will follow the procedure described in [TH4T6] . Moreover, the discussion is based on 
the assumption that the effective particle-wall interaction time is small compared to any 
characteristic time of interest and no permanent adsorption occurs jH]. 

For the high-order models, as "particles" from more than one layer of computational grids 
can hit the wall, we have to properly identify them in order to implement the boundary 
condition. For this purpose, N layers of ghost grids are introduced (see the example of the 
D2Q17 lattice and its grid arrangement shown in Fig|l]and[2]), where N can be determined 
via the corresponding maximum value of the discrete velocity heading towards the wall (e.g., 
N = 3 for the D2Q17 lattice). As a common practice, the physical wall is located at the half 
grid space between the ghost and fluid grids. To further distinguish incoming and outgoing 
particles, we use d a x and c a j to represent their velocities respectively, where / denotes the 
layer number of the ghost grid ranging from to AT — 1. Similarly, the distributions of 
incoming and outgoing particles at layer I are written as fai(x w ,t) and /^(cc^t), where 
the superscripts I and O stand for 'incoming' and 'outgoing'. The corresponding discrete 
velocities must satisfy the condition (c' a t — u w ) ■ ndt < —Idx and (c a j — u w ) ■ ndt > Idx, n 
denotes the unit vector normal to the wall surface dfl at x and directed from the wall into 
the gas. Note in the present lattice system dx/dt = 1 so the conditions are equivalent to 
(c' a l — u w ) n < —I and (c a j, — u w ) n > I. Indeed, c' a l and symmetric pair. On the 

other hand, the known information of the wall, i.e., the position, velocity and temperature, 
are represented by x w ,u w and T w . 

Obviously, the distribution f^ t (x w ,t) can be obtained by naturally streaming the distri- 
bution function at fluid grids into the corresponding ghost ones. We need to determine the 
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unknown distribution f^ d (x w ,t) according to the principle of diffusion reflection. Similar to 
the derivation of the continuum version of diffusion-reflection condition |14j . we first write 
down the mass of outgoing and incoming particles as, 



M°i = f^(x w ,t)dV xEdQ,(c a!l -u w )-n>l, (6) 



M T a j = fij{x w ,t)dV x G dQ, (c' al — u w ) ■ n < —I, (7) 

where Ai stands for mass and dV denotes the volume of the grid cell. It is worth noting 
again here that, due to the exact advection of the LB method (cf. Eq.(|5])), the flux term in 
the continuum version (cf. Eq. (1.11.1) in [14J), can be replaced by the distribution function 
itself. Hence, according to the mass conservation, we have, 
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where R(c' al — >■ c a j,x w ,t) is the so-called scattering probability. Immediately, we arrive at 
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Moreover, the scattering probability R must satisfy the property of non-negativeness, nor- 
malization and reciprocity condition[H]. Particularly, the normalization condition, corre- 
sponding to mass conservation under the assumption of no permanent adsorption, can be 
written as, 

N-1 

E R{d al ^c ajl ,x w ,t) = l. (10) 

(=0 (c a j—u w )-n>l 

So far, the discussion is still generic as we have not introduced any specific assumption for 
the diffuse- reflection principle. Therefore, the above formulation may also be used to derive 
other type of boundary condition. 

If the assumption of the diffuse-reflection boundary condition is applied, the scattering 
probability can be easily calculated as 



S/Iq 1 Y.(c a ^-u w )-n>l9a(u w ,T u 
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Figure 1: Illustration of the D2Q17 lattice. Each discrete velocity is represented by the length and 
direction of the line connecting the origin of coordinate and the corresponding dotted point. 



Hence, the distributions of outgoing particles can be written as 
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IV. NUMERICAL VALIDATION 

To validate the proposed implementation of kinetic boundary condition, we consider 
steady Couette flow confined in two parallel planar plates located at Y = and Y = 1 
and moving oppositely with the same speed. All the quantities are presented in their non- 
dimensional form, and both isothermal and thermal conditions are considered. Therefore, 
three lattice systems, namely D2Q17[1], D2Q16[l7] and D3Q121[12l[T8] are tested, where the 
D2Q17 and D2Q16 models are appropriate for the isothermal cases and the D3Q121 model 
for the thermal ones. The D2Q17 lattice is illustrated in FigjTJand the corresponding grid 
arrangement is shown in Fig|2j The details of three lattices are omitted here for simplicity, 
which can be found in pQ, [T7] and [Tg] . 

We restrict to the Couette flow within the slip-flow regime, so we may be able to use 
solutions of the Navier-Stokes-Fourier (NS) equations as reference. For the NS solutions, it 
is necessary to apply the velocity-slip and temperature-jump boundary conditions so that 
the velocity and temperature profiles can be written as 
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Figure 2: Illustration of the arrangement of the ghost grid and the corresponding outgoing discrete 
velocity (arrow at bottom for the D2Q17 model (see Fig{TJ. The point denotes the bulk grid while 
the dashed circle represents the ghost grid. 

and 

Kn (8C P T W + hul) + 8C p T w Kn 2 + 2C P T W - APrU w (Y - l)Y 

2C P (2Kn + 1) 

where the Knudsen number is denned as 

Kn = , /J^M. (15) 

U w denotes magnitude of the component of interest of the wall velocity u w . For some 
relatively larger Knudsen numbers we may also compare to the solution of the linearized 
Boltzmann-BGK (L-BGK) equation. 

We first evaluate the D2Q17 and D2Q16 models for isothermal flows which are presented 
in Fig|3]and[4j where U w is set to be 0.05. Both models are simulated with 100 computational 
grids in the direction of interest and the comparisons are made against the NS solutions for 
Kn < 0.05 and the L-BGK solutions for Kn > 0.05 respectively. The results show that 



the boundary condition Eq.(12) works correctly for the isothermal flows. For Kn < 0.05, 
the velocity profiles are captured well by the D2Q17 model while some deviations from the 
L-BGK results are observed for larger Knudsen numbers, particularly at Kn = 0.1 (see 
Fig|3]). However, this is of no surprise as it is known that these deviations are due to the 
lattice structure pEj. A further comparison to the finite difference (FD) implementation 
of Eq.Q (see the description in [19], where the numerical simulation is validated by the 
analytical solution in [20]), confirms the appropriateness of the boundary implementation. 
Interestingly, the D2Q16 model can given much better predictions for the velocity profile. 
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Figure 3: The velocity profiles for the isothermal cases with the D2Q17 model. The velocity is 
further normalized by the wall velocity. 
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Figure 4: The velocity profiles for the isothermal cases with the D2Q16 model. The velocity is 
further normalized by the wall velocity. 

Even at Kn = 0.5, it still gives satisfactory results, see Fig. |4j The reason was already 
discussed in [T9"] 

For further validation, we also simulate the thermal Couette flows using the D3Q121 
model. For these flows the relevant parameters are C p = 5/2 and Pr = 1 while the wall 
temperatures are set to be 1 and their speed U w is set to be 0.2. As relatively small Knudsen 
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Figure 5: The temperature profiles for the thermal cases with U w = 0.2. 

numbers are considered here, the results are compared to the NS solutions, see Figj5j The 
subtle temperature jumps are well captured at the wall boundary. These agreements again 
confirm the appropriateness of the boundary treatment. The velocity profiles show similar 
behavior to the isothermal cases, so they are not presented in Fig. [5j 

To evaluate the numerical accuracy, a convergence study is conducted for the ther- 
mal case of Kn = 0.01. The simulations are run for six different grid resolutions 
N G = 20,50,100,200,500,1000 in the direction of interest. The results of N G = 1000 is 
then chosen as reference and the global relative errors of the velocity and temperature are 
defined as 

nil = \ _/V„ ,„s > l lb J 



and 



En 



where Ua and Ta represent the results of N G 
accuracy is achieved globally 



EgiK-^O-rA^)] 2 (17) 

1000. FigM shows that the second order 



V. CONCLUDING REMARKS 



To conclude, we have formulated the kinetic diffuse reflection boundary condition for 
high-order LB models with emphasis on retaining the "streaming-collision" mechanism. The 
numerical tests for both isothermal and thermal Couette flows show that the present bound- 
ary condition can capture velocity-slip and temperature-jump very well within the capacity 
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Figure 6: Dependence of the Ejj and Et on the grid number Nq- 

of the corresponding lattices. In term of numerical accuracy, we show that the second order 
accuracy can be achieved globally. 
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